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We study the Hubbard-Holstein model, which includes both the electron-electron and electron- 
phonon interactions characterized by U and g, respectively, employing the dynamical mean-field 
theory combined with Wilson's numerical renormalization group technique. A zero temperature 
phase diagram of metal-insulator and insulator-insulator transitions at half-filling is mapped out 
which exhibits the interplay between U and g. As U (g) is increased, a metal to Mott-Hubbard 
insulator (bipolaron insulator) transition occurs, and the two insulating states are distinct and can 
not be adiabatically connected. The nature of and transitions between the three states are discussed. 



The interaction-induced phenomena are a fundamental 
problem because of the new physics that emerges. For in- 
stance, the interaction-driven metal-insulator transition 
(MIT) has been studied to a great extent for systems 
which involve the electron-electron or electron-phonon 
interactions [ij. In real condensed-matter systems, both 
interactions exist and it will be important to understand 
the interplay between them. The prototype model for 
describing such interplay is the Hubbard-Holstein (HH) 
model, which includes the onsite electron-electron and 
electron-phonon interactions characterized by U and g, 
respectively, as given by Eq. |^ below. The Hubbard 
model {g — 0) and Holstein model ([/ = 0), which are 
special limits of the more general HH model, have been 
extensively studied using the dynamical mean-field the- 
ory (DMFT) in infinite dimensions in the context of 
the interaction-driven MIT. Both models exhibit MIT as 
[7 or 5 is increased above a critical value. Now, consider 
the MIT in the U — g plane. Natural questions then arise 
as to whether and how the two insulating states, metal- 
lic states, and metal-insulator transitions of the Hub- 
bard and Holstein models are different, and how they 
are affected when both U and g are present and compete 
with each other. These questions will be addressed here 
by studying the HH model within the DMFT combined 
with Wilson's numerical renormalization group (NRG) 
01 . The NRG technique 5, 6] is particularly powerful 
in that it is non-perturbative in nature so that it can 
cover the whole parameter space and that it can probe 
the extremely small energy scales like the narrow coher- 
ence peak and soft phonon mode which emerge as the 
metal-insulator transitions are approached. 

We will present a zero temperature phase diagram of 
unbroken symmetry ground states with a focus on the 
nature of and transitions between the ground states as U 
and g are varied. The main results of the present work 



are: (i) The ground state is a metallic (M) state when 
both U and g are small, but is a bipolaron insulating 
(BPI) state when g is large, and is a Mott-Hubbard in- 
sulating (MHI) state when U is large (see Fig. 1). In 
the shaded region in Fig. 1, the insulating and metallic 
solutions coexist, and the ground state is the one with a 
lower energy, (ii) The phase transition between M and 
BPI occurs along the dashed line that cuts through the 
coexistence region connecting E and T of Fig. 1. The 
order parameter {{n — l), see below.) changes discontin- 
uously at the transition and the M-BPI transition is first 
order (Fig. 2). (iii) The transition between M and MHI 
occurs along the outer boundary of the coexistence re- 
gion denoted by the solid line connecting T and U2. It is 
second order as in the Hubbard model (Fig. 3). (iv) The 
two kinds of insulators (BPI and MHI) are distinct and 
can not be adiabatically connected. The transition be- 
tween them is first order, and occurs along U — 2g^/wo 
(Fig. 4. See below). These points will be discussed in 
more detail below. 

The HH model is defined by 

^ <i,j><J i 

-\- luq^ ajfli + ^ (a] + flj^ (rii - I) , (I) 

i i 

where q is the nearest-neighbor coordination number, 
n-ia — cl^Cia the electron density operator for spin a at 
site i, and rii = ^^riia-- The electrons have the onsite 
Coulomb repulsion U, and are linearly coupled with the 
Einstein phonon of frequency wq with the onsite coupling 
constant g. The DMFT was employed to solve the HH 
model of Eq. ^ . The DMFT maps a lattice model onto 
an effective single site impurity model imbedded in a bath 
which is determined self-consistently via iterations 0, Q . 
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FIG. 1: The zero-temperature metal-insulator transition 
phase diagram of the HH model at half-filling. The cuts la- 
beled with 2 (a & b), 3, and 4 correspond, respectively, to 
the Figs. 2, 3, and 4. The cuts 2a, 2b, 3, and 4 represent, re- 
spectively, the calculations along the fixed g = 0.16, U = 1.0, 
g = 0.1, and U = 2.5. The shaded area is the coexistence 
region where both the metallic and insulating solutions ex- 
ist. The dashed lines between M and BPI, and between BPI 
and MHI states represent 1st order phase transitions, and the 
solid line between M and MHI represents 2nd order. See the 
text for details. 



The effective impurity problem was solved by the NRG 
technique in the present study. 

The NRG-fDMFT has been successfully appHed to 
the half-filled Hubbard model to study the MIT at both 
T = 0] and finite temperatures , and to the half- filled 
Holstein model at T = . For the Hubbard model, the 
ground state is insulating for U > Uc2 and metallic for 
U < Uci- Between Ud < U < Uc2, the metallic and 
insulating solutions coexist, and the metallic state has 
the lower energy and is the ground state. The MIT oc- 
curs at Uc = Uc2, and is 2nd order. For the Holstein 
model, the ground state is insulating for g > gc2 and 
metallic for g < gd ■ Meyer et al. reported that gc2 — gci 
is reduced as the phonon frequency ujq is decreased and 
gel = gr2 for (jjQ = 0.05 W, where W = At is the band- 
width y. We generalized these lines of research and 
study the HH model at half-filling with the NRG-f DMFT 
incorporating the improved method for calculating elec- 
tron and phonon spectral function 10]. We adopt the 
semi-circular density of states which is realized in the 
Bethe lattice in infinite dimensions. 

Our central result is the Fig. 1, where we present a 
zero temperature phase diagram in the U — E^p plane 
of unbroken symmetry ground states of the HH model 
at half-filling obtained from the NRG-^DMFT calcula- 
tions. Ei,p = 2g'^ /loq is the bipolaron binding energy 
which is twice the polaron energy, E^p = 2Ep. We take 
the bandwidth W as the unit of energy {W = 1). We 



take ujQ = 0.05 as with Meyer et al. Three distinct states 
are found in the U — E^p plane: metal, bipolaron insula- 
tor, and Mott-Hubbard insulator states. We will look at 
them in more detail below. Specifically, we will address 
(i) the nature of M, BPI, and MHI states, (ii) transi- 
tion between M and BPI, (iii) transition between M and 
MHI, and (iv) transition between the BPI and MHI. The 
results obtained from the NRG-I-DMFT regarding these 
points were summarized in the introduction above. 

(i) Nature of M, BPI, and MHI states. - The ground 
state is a metallic state when both U and g are small. 
The M state becomes more correlated as g or U is in- 
creased as reflected in the decreasing quasiparticle weight 
z = [1 — d'Ei{uj)/duj]~l^f^ , where the subscript 1 (2) 
refers to the real (imaginary) part, z vanishes at the 
outer boundary of the coexistence region. The BPI state 
is insulating because the attractive interaction between 
the electrons mediated by the phonons binds two elec- 
trons into a bipolaron as Capone and Ciuchi reported 
The MHI state, on the other hand, is insulating 
because of the repulsion U between two electrons at the 
same site. 

An important point for understanding the overall fea- 
ture of the phase diagram is that the effective interaction 
between two electrons in the HH model is given, after in- 
tegrating the phonons out, by 
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(2) 



The interactions between electrons due to g and U com- 
pete each other for w ^ wq, which yields Veff{uj = 0) = 
along U = Ehp. The U = Ehp line plays a special role 
in the phase diagram in that all three states meet at the 
triple point T lying on this line and the BPI-MHI tran- 
sition occurs along it, as we will discuss below. Also, the 
critical gc, which separates M and BPI along E and T, 
should increase as U is increased because the attractive 
interaction due to g has to overcome the repulsive inter- 
action of U to bind two electrons into a bipolaron. This 
expectation is indeed borne out as shown in Fig. 1. On 
the other hand, the Uc, which is the solid line connecting 
T and U2, depends weakly on g. 

For the coexistence region with U > Ei,p, the argu- 
ments which led to the conclusion of the MIT of the 
Hubbard model being a 2nd order at U = Uc2 still 
apply and the M is the ground state. On the other hand, 
in the coexistence region with E^p > J7, the {i^ii^i) gpj 
and (n|nx)j^^, where {O) j^^^j^ stands for an expectation 
value of an operator O in the metallic (insulating) state, 
do not meet as shown in Fig. 2(a). This implies that the 
M-BPI transition is 1st order. The energy calculations in 
Fig. 2(a) show that the M and BPI solutions cross each 
other inside the coexistence region. 

(ii) Transition between M and BPI. - The MIT of half- 
filled Holstein model in infinite dimensions was discussed 
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FIG. 2: (a) The ((n - l)^)/2, which is equal to the double 

occupancy (niUiJ in M state, and ground state energy Eg are 
plotted along the cut 2a with a fixed g = 0.16 in Fig. 1 as U 
is increased (for M) or decreased (for BPI). It can be read off 
from the plots that the coexistence region is 0.925 < U < 1.0, 
and the phase transition occurs around U — 0.96 where the 
energies of the M and BPI states become the same, (b) Plots 
of z, C, and ^(a) as g is increased along the cut 2b of Fig. 
1. In Figs, (c) and (d), the electron and phonon spectral 
functions, A{u!) and B{uj), respectively, are shown along the 
same cut with (b). 



previously by Meyer et al. [3 and by Benedetti and Zey- 
her 13 . The BPI state is insulating because electrons are 
bound in bipolarons. Because we set the chemical poten- 
tial /X = U/2, which prefers the electron density (n) = 1, 
the bipolaron formation may be accommodated by re- 
constructuring the system into a phase separated state 
or a charge ordered state where the doubly occupied and 
empty sites alternate in the real space. This possibil- 
ity is not allowed in the present DMFT work because 
the unit cell consists of a single site. Instead, the bipo- 
laron instability shows up as degenerate ground states: 
In the NRG-I-DMFT calculations, the ground state is 8 
fold degenerate in the BPI regime, with a one-to-one cor- 
respondence between the two sets of (n) — 1 w 1 and 
(n) — 1 w —1, which on averaging satisfies (n) = 1. The 
order parameters for the M-BPI may be taken as 



C={n-l). 



(3) 



Now, consider the cut 2b in Fig. 1 along which we in- 
crease g with U = 1.0. Let gd and gc2 be, respectively, 
the lower and upper boundaries of the shaded coexistence 
region. As g — > gc2, the quasiparticle weight z and the 
renormalized phonon frequency Q approach continu- 
ously, but, the order parameter C increases from to a 
finite value discontinuously, as shown in Fig. 2(b) and 
(d). The softening phonon mode is a manifestation of 
a lattice instability just like structural phase transitions. 
A stability is restored by a condensation of the unstable 
mode. It results in a non-zero expectation value of the 
phonon operator {{^aj ^ 0), which may be taken as an 
order parameter for the M-BPI transition. The phonon 
hardens back to the bare mode as g is increased above 
gel because the screening is not effective in an insulating 
state. The a and (n—X) are linearly coupled in the Hamil- 
tonian of Eq. J^l, and a zero/nonzero (a) is expected to 
lead to a zero/nonzero C. Therefore, an equally good 
order parameter is the C of Eq. (j^J. 

In Fig. 2(b), we show z, C, and ^(a) with increasing 
g. Note that (n — l) « ^('^) ^ expected from a simple 
mean- field theory. In Fig. 2(c), we plot the renormal- 
ized electron spectral function, A(w) = —^Im(lc,c^))^, 
to show the MIT as g is increased along the cut 2b. In 
Fig. 2(d), we plot the renormalized phonon spectral func- 
tion, B{uj) — ~^Im{{a+a\ a+a^))uj, calculated as g is in- 
creased using the method in Note that a.t g = 0.1664 
which is close to gc2 = 0.1665, the phonon mode splits 
into two components as reported by Jeon et al. for the 
Anderson-Holstein model: one component develops into 
the soft mode and the other hardens back to the bare 
mode But, in the coexistence region where gd and 
gc2 are not close, the 1st order MIT from M to BPI pre- 
empts the emergence of the soft mode as g is increased 
because the MIT occurs at g^ (g^i < .9c < 9c2)- 

The M-BPI transition may be understood from the 
Ginzburg-Landau theory. The free energy may be writ- 
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FIG. 3: In Fig. (a), we plot as C/ is increased, and (n-\ni^ 
as U is increased and decreased along the cut 3 of Fig. 1. In 
Fig. (b), the electron spectral functions Aiyi) are shown as U 
is increased along the same cut. 

ten in terms of the the local phonon coordinate ^ = 

F +\ie , (4) 

with 7 > 0. The phase transition can be 1st order for /? < 
0. In the context of the M-BPI transition, the gc\ < 9 < 
gc2 regime corresponds to /3^/47 > a > 0, for which three 
local minima exit as a function of ^. The gd, gc, and gc2 
correspond, respectively, to a = /3^/47, 3/3^/167, and 0. 
Although the a and f3 are complicated functions of U and 
g, the DMFT phase diagram of Fig. 1 indicates that a 
and P are mainly determined by, respectively, g and U. 
This observation suggests an interesting possibility that 
the M-BPI transition may become 2nd order for small 
values of U where P > 0, and a quantum tricritical point 
exists between E and T. This possibility will be explored 
further in a subsequent study. 

(in) Transition between M and MHI. - As the MIT 
from M to BPI is driven by growing charge suscepti- 
bility, the one from M to MHI is driven by growing 
spin susceptibility. In Fig. 3(a), we plot z as J7 is in- 
creased, and (ri|n|) as U is increased and decreased 
for a fixed g = 0.1. z goes to continuously, and 
{n-T^ni) — (iT'^iT'l) MHI as C/ — > Uc2- The M-MHI tran- 



FIG. 4: The order parameters C is plotted as a function of 
g for the fixed U = 2.5 corresponding to the cut 4 of Fig. 
1. The phase transition between BPI and MHI occurs along 
U = Ebp. 

sition at Uc = Uc2 is, therefore, 2nd order. The phonons 
are hardly renormalized (not shown) because the charge 
degree of freedom is not soft unlike the MIT between M 
and BPI. In Fig. 3(b), A{uj) is plotted to show the MIT 
as U is increased. 

The Uc does depend on g, although weakly, as can 
be seen from the Fig. 1, which implies that the M-MHI 
transition is coupled with the phonons. This raises an in- 
teresting possibility that some discrepancies between the 
experimental observations on V2O3 systems and DMFT 
results of the Hubbard model mav be resolved in terms 
of the electron- phonon coupling |l4| . 

(iv) Transition between BPI and MHI. - Now, we turn 
to an interesting question about the two kinds of insu- 
lators: Are BPI and MHI different, or, put differently, 
can the two insulating states be adiabatically connected? 
Along the line U = E^p, the effective interaction in the 
static limit vanishes. One might then naively expect that 
a metallic state is the ground state along this line all the 
way to U — Ehp 00. One can show, however, that for 
U — Ehp 00 the ground state is an insulating state 
using the exact atomic limit Green's function [lolll5l |. as 
found from the detailed NGR-I-DMFT calculations. 

In Fig. 4, we show C as 5 is varied for U — 2.5. The 
C changes discontinuously exactly at E^p = U and does 
not exhibit any other non-analytic behavior as g is varied. 
We checked another point on the line E^p — U by varying 
U with a fixed g = 0.3. The transition occurs exactly 
at E\,p = U (not shown). This confirms that the two 
insulating states are distinct, and the phase transition 
between them occurs along the line of C/ = E\yp and is 
1st order. The point T in Fig. 1 is a triple point where 
two 1st order and one 2nd order phase transition lines 
meet. 
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